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Abstract 

The  next-generation  all-electric  ship  represents  a  class  of  design  and  control  problems 
in  which  the  system  is  too  large  to  approach  analytically,  and  even  with  many  conven¬ 
tional  computational  techniques.  Additionally,  numerous  environmental  interactions 
and  inaccurate  system  model  information  make  uncertainty  a  necessary  consideration. 
Characterizing  systems  under  uncertainty  is  essentially  a  problem  of  representing  the 
system  as  a  function  over  a  random  space.  This  can  be  accomplished  by  sampling 
the  function,  where  in  the  case  of  the  electric  ship  a  ”  sample”  is  a  simulation  with 
uncertain  parameters  set  according  to  the  location  of  the  sample.  For  systems  on  the 
scale  of  the  electric  ship,  simulation  is  expensive,  so  we  seek  an  accurate  representa¬ 
tion  of  the  system  from  a  minimal  number  of  simulations.  To  this  end,  collocation 
is  employed  to  compute  statistical  moments,  from  which  sensitivity  can  be  inferred, 
and  to  construct  surrogate  models  with  which  interpolation  can  be  used  to  propagate 
PDF’s.  These  techniques  are  applied  to  three  large-scale  electric  ship  models.  The 
conventional  formulation  for  the  sparse  grid,  a  collocation  algorithm,  is  modified  to 
yield  improved  performance.  Theoretical  bounds  and  computational  examples  are 
given  to  support  the  modification.  A  dimension-adaptive  collocation  algorithm  is  im¬ 
plemented  in  an  unscented  Kalman  filter,  and  improvement  over  extended  Kalman 
and  unscented  filters  is  seen  in  two  examples. 
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Chapter  1 


Introduction 


Large-scale  engineered  systems  pose  a  class  of  design  problems  for  which  there  is 
currently  no  existing  comprehensive  approach.  The  space  of  possible  realizations  can 
be  immense  for  such  problems  due  to  dimensionality,  placing  them  well  outside  the 
computational  limits  of  traditional  techniques.  The  all-electric  ship  is  a  representa¬ 
tive  example  of  large-scale  systems,  and  its  development  is  of  significant  interest  to 
the  control  and  simulation  communities  [9,23,35].  Its  superiority  over  current  ships 
includes  high  efficiency  power  transfer  through  the  ship,  rapid  reconfiguration,  layout 
flexibility,  and  azimuthing  thrusters  providing  enhanced  maneuvering  and  hydrody¬ 
namic  efficiency;  presently,  this  is  only  a  potential.  To  find  designs  and  corresponding 
control  strategies  that  attain  this  potential,  the  system  must  be  considered  as  a  whole; 
from  this  perspective,  it  is  the  size  of  the  problem  that  defines  it. 

Very  few  real  systems  can  be  truly  described  by  a  mathematical  model.  Reasons 
for  this  include  system  parameters  that  can  only  be  known  to  a  certain  accuracy, 
unpredictable  interaction  with  the  environment,  and  dynamics  unaccounted  for  in 
model  equations.  Large-scale  systems  have  numerous  components  and  many  modes 
of  interaction  with  the  environment,  and  consequently  tend  to  have  many  sources  of 
uncertainty. 

A  standard  way  to  broaden  the  scope  of  a  model  is  to  represent  uncertain  aspects 
of  the  actual  system  with  random  variables  or  processes,  resulting  in  systems  of 
stochastic  differential  equations  [22,29].  Each  random  variable  in  a  model  contributes 
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a  dimension  to  the  sample  space  of  the  model,  which  is  the  space  of  all  possible 
realizations.  Large  systems  often  have  high-dimensional  sample  spaces,  and  because 
they  are  usually  analytically  intractable,  they  suffer  what  is  known  as  ’The  Curse  of 
Dimensionality’,  a  term  coined  by  Richard  Bellman  to  describe  the  computationally 
difficulties  that  arise  in  high-dimensional  spaces. 

In  this  work,  collocation  algorithms  for  numerically  solving  stochastic  differential 
equations  are  explored.  Information  is  extracted  from  the  resulting  solutions  in  the 
form  of  statistical  moments  and  probability  distribution  functions  (PDF’s),  which  are 
respectively  computed  via  numerical  integration  and  interpolation.  These  tools  are 
applied  to  electric  ship  models  and  unscented  Kalman  filtering. 

This  thesis  makes  three  basic  contributions.  First,  a  framework  for  analyzing 
power  system  models  under  uncertainty  is  developed  and  applied  to  three  different 
systems,  a  notional  electric  ship  model,  a  pulse  power  system,  and  a  larger  notional 
electric  ship  model  undergoing  a  charging  event.  For  each  system,  sensitivity  is  in¬ 
ferred  from  a  different  perspective.  Second,  a  modification  to  the  the  sparse  grid,  a 
collocation  algorithm,  is  presented  and  justified  via  theoretical  bounds  and  numerical 
experiments.  Third,  a  dimension-adaptive  collocation  algorithm  is  then  implemented 
within  the  unscented  Kalman  filter,  and  improvement  is  found  over  other  Liters  on 
two  example  systems. 


16 


Chapter  2 


Uncertainty  in  Dynamic  Systems: 
Setting  the  Stage 


Uncertainty  can  be  incorporated  into  dynamical  systems  such  as  the  electric  ship 
by  modeling  them  with  stochastic  differential  equations.  In  this  chapter,  a  cursory 
overview  on  stochastic  differential  equations  and  methods  for  solving  them  numer¬ 
ically  is  given,  followed  by  three  examples  of  dynamical  systems  with  uncertainty, 
the  latter  of  which  is  the  all-electric  ship.  From  a  general  perspective,  we  seek  the 
solution  x  :  fl  x  9?+  — »  K  to  the  stochastic  differential  equation 

L(t,uj;x)  =  f(t,Lu),  Xo  =  g(uj),  1  6  9f+,  c v  E  fi,  (2.1) 

where  L  is  the  differential  operator  and  the  sample  space  [45].  Equation  (2.1)  may 
contain  both  random  variables  and  random  processes.  A  random  process  may  be 
represented  by  a  random  variable  at  each  moment  in  time;  this  however  makes  the 
sample  space  U  infinite  dimensional,  a  feature  that  cannot  be  directly  accommodated 
computationally.  The  Karhunen-Loeve  expansion  [14]  offers  a  means  to  approximate 
random  processes  using  a  finite  number  of  random  variables  in  a  truncated  spec¬ 
tral  expansion.  Unfortunately,  accurately  representing  a  random  process  tends  to 
dramatically  increase  the  dimension  of  U. 

For  all  but  the  simplest  of  systems,  (2.1)  must  be  solved  numerically.  There 
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are  multiple  approaches,  including  brute  force  Monte  Carlo  simulation,  polynomial 
chaos  [42],  and  collocation,  the  latter  of  which  we  focus  on  here. 

Useful  information  is  generally  extracted  from  a  solution  x(t,ou)  to  (2.1)  through 
integration,  i.e.  finding  mean  and  variance  trajectories,  and  interpolation,  which 
enables  the  construction  of  PDF’s. 

The  expected  value  of  a  continuous  random  variable  c o  and  that  of  a  function  of 
a  continuous  random  variable  g( u)  are  given  respectively  by 

E[ uj\  =  [  u>p(u>)duj,  E\g(uj)\  =  f  g(co)p(co)du},  (2.2) 

Jn  Jn 

If  we  view  equations  in  the  form  of  (2.1)  as  functions  of  random  variables,  moments 
of  the  trajectories  can  than  be  calculated  as  in  (2.2): 

E[x(t,u>)]  =  [  x(t,u)p(u)du.  (2.3) 

Jn 

x(t,u)  is  the  solution  to  an  ODE  corresponding  to  a  specific  realization  in  O.  Rarely, 
x(t,u)  may  be  calculated  analytically;  more  often  it  is  obtained  through  simulation. 
For  large-scale  systems,  simulations  can  be  computationally  very  expensive,  so  it  is 
desirable  to  employ  an  efficient  discretization  of  (2.3),  i.e.  one  in  which  high  accuracy 
is  achieved  from  a  minimal  number  of  simulations.  In  the  next  three  sections  examples 
of  dynamical  systems  with  uncertainty  are  given. 

The  examples  in  this  chapter  illustrate  how  uncertainty  is  incorporated  into  dy¬ 
namical  systems  and  the  information  that  is  gained  from  their  solutions. 


2.1  Uncertainty  in  a  First  Order  System 

This  example  demonstrates  simply  the  calculations  carried  out  on  more  complicated 
systems  later  in  this  thesis.  Consider  the  first  order  stochastic  ordinary  differential 
equation 

x  +  RiX  —  0,  x(0)  —  R2,  (2-4) 
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where  R\  and  R2  are  uniformly  distributed  random  variables  in  the  interval  [0  1]. 
The  analytical  solution  is  given  by 

x(t,  R)  =  R2e~Rlt.  (2.5) 

The  mean  trajectory  can  then  be  computed  analytically  by  integrating  (2.5)  multi¬ 
plied  by  one  (because  R  is  uniform)  over  the  system’s  random  space,  [0  1]  x  [0  1]: 

E[:r!/,Rf  =  [  [  r2e~Tltdridr2  —  j-(l  —  e-i). 

Jo  Jo  2t 

The  variance  can  be  computed  similarly,  and  is 

var(x(t,  R))  —  [  [  (roe^rit  —  E[x(t)])2  dr\dr2  =  — (1  —  e~2t) 

Jo  Jo  6 1 

Notice  that  because  this  system  has  two  random  parameters,  the  moment  calculations 
are  double  integrals. 


-uE-e-V- 


2.2  Uncertainty  in  a  Double  Pendulum 

Now  consider  a  double  pendulum  with  uncertainty  in  the  initial  displacement  of  each 
arm  (Fig.  2-1).  As  in  the  previous  example,  this  system  has  a  two  dimensional 
random  space.  However,  the  double  pendulum  equations  (see  [37])  cannot  be  solved 
analytically,  and  so  (2.3)  must  be  solved  numerically  for  this  system,  both  in  the  sense 
that  the  solution  of  the  system  equations  x(f ,  R)  must  be  obtained  through  simulation 
and  that  the  statistical  moment  integrals  must  also  be  computed  by  discretizing  the 
random  space. 

Fig.  2-2  shows  the  mean  trajectories  with  standard  deviation  envelopes  for  this 
system.  The  solutions  to  the  system  equations  were  generated  using  the  fourth  order 
Runge-Kutta  scheme,  and  the  moments  were  then  computed  using  techniques  de¬ 
scribed  in  chapter  3.  Notice  that  at  about  2.3  seconds  the  variance  grows  very  large. 
This  is  because  the  system  after  that  time  is  sensitive  to  the  prescribed  uncertainty, 
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Figure  2-1:  Double  pendulum  with  uniform  uncertainty  in  the  initial  positions. 

meaning  that  it  can  exhibit  a  wide  range  of  behavior  depending  on  where  in  the  ran¬ 
dom  space  the  system  is  realized.  Mathematical  definitions  of  sensitivity  will  be  given 
in  chapter  5,  and  will  be  used  to  quantify  the  robustness  of  systems  to  uncertainty. 

2.3  Uncertainty  in  the  All-Electric  Ship 

Physically,  uncertainty  in  the  all-electric  ship  comes  in  many  forms.  Values  of  physical 
parameters  like  capacitances,  inertia  constants,  and  flux  leakages  are  only  known  to 
manufacturers’  specifications,  and  may  change  with  time  and  use.  There  is  significant 
environmental  interaction  in  the  form  of  wind  and  waves;  a  large  wave  can  partially 
of  fully  expose  the  propellers,  dramatically  and  suddenly  altering  the  torque  load 
on  the  induction  motors.  Lastly,  there  is  a  diverse  range  of  operating  conditions, 
with  each  mode  utilizing  the  ship  systems  differently.  Traveling  long  distances  will 
place  the  ship  under  fairly  steady  conditions,  while  in  a  combat  scenario  there  will 
be  potentially  heavy  loading  from  propulsion  and  radar,  irregular  pulsed  loads  from 
weapon  systems,  and  the  ultimate  uncertainty,  damage.  Clearly,  modeling  the  full 
scope  of  uncertainty  that  might  be  seen  by  the  all-electric  ship  (or  any  large-scale 
system)  would  take  one  well  outside  of  the  framework  presented  here. 

Simplified  analysis  is  revealing  nonetheless,  and  often  simplifying  one  aspect  of 
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Mean  solutions  with  standard  deviation 


Figure  2-2:  Mean  solutions  with  standard  deviation  envelopes  for  a  double  pendulum 
with  uniform  initial  position  uncertainty. 


an  analysis  allows  greater  depth  in  another.  By  restricting  the  types  of  uncertainty 
considered,  the  size  of  the  system  one  can  study  is  increased.  Plainly,  by  allowing 
only  parametric  and  initial  condition  uncertainty  in  models  of  large-scale  systems, 
it  becomes  possible  to  study  the  effect  of  that  uncertainty  not  on  components  or 
portions  of  the  systems,  but  on  the  systems  as  a  whole.  In  chaper  5,  such  analysis  is 
applied  to  the  all-electric  ship. 
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2.4  Summary 


In  this  chapter,  a  brief  overview  of  how  uncertainty  i: 
modeled  with  stochastic  differential  equations  is  given, 
were  shown. 


dynamical  systems  can  be 
Three  illustrative  examples 
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Chapter  3 


Technical  Overview 


In  this  chapter,  collocation  algorithms  for  solving  stochastic  differential  equations 
are  presented,  along  with  basic  elements  from  which  they  can  be  constructed  and 
methods  for  extracting  information  from  the  resulting  solutions. 

Collocation  approximates  solutions  to  (2.1)  by  calculating  solutions  to  determin¬ 
istic  differential  equations  corresponding  to  special  points  in  the  random  space 
and  then  using  those  solutions  to  construct  a  function  over  If  at  each  moment  in 
time.  Monte  Carlo  methods  differ  only  in  the  choice  of  points,  using  pseudo-  or 
quasi-random  points  rather  than  collocation  points. 

At  this  point,  it  is  sensible  to  simplify  the  discussion  from  stochastic  differential 
equations  to  functions;  a  solution  to  a  stochastic  differential  equation  at  a  fixed  time 
is  essentially  a  function  over  a  random  space,  and  in  the  case  of  collocation,  finding 
that  solution  amounts  to  approximating  a  function  at  each  moment  in  time.  Two 
basic  settings  in  which  collocation  algorithms  are  used  are  numerical  integration  and 
interpolation  of  functions. 

3.1  Numerical  Integration 

Probability  says  that  means,  variances,  and  higher  moments  of  random  variables 
are  integrals.  These  integrals  must  be  computed  numerically  for  large-scale  systems, 
and  thus,  because  of  ’The  Curse  of  Dimensionality’,  become  problems  of  efficiently 
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sampling  a  multi-dimensional  random  space;  this  is  a  quintessential  aspect  of  this 
work. 

Numerical  integration  [6]  is  essentially  a  problem  of  efficiently  characterizing  a 
function  over  a  space  by  sampling  it  at  discrete  points.  Considered  here  are  algorithms 
that  can  be  expressed  as  a  sum  of  function  evaluations  multiplied  by  weights,  such 
that 

J /(x)dx  ~  Y,  wif  (3.1) 

i 

Later,  these  algorithms  are  applied  to  (2.3)  to  compute  moments  of  trajectories  of 
dynamical  systems. 

3.1.1  In  One  Dimension 

Familiar  one-dimensional  numerical  integration  schemes  (also  called  rules)  are  Ric- 
mann  sums,  the  trapezoidal  rule,  and  Simpson’s  rule.  These  methods  produce  rel¬ 
atively  crude  results  in  comparison  with  polynomial-based  approaches,  but  have  in 
common  that  they  can  all  be  written  in  the  form  of  (3.1),  and  thus,  once  the  weights 
and  function  evaluation  points  (sometimes  called  nodes)  have  been  computed,  are 
identical  in  implementation.  The  set  of  points  and  weights  comprising  the  level  i 
rule  of  a  numerical  integration  scheme  will  be  denoted  $*.  With  this  notation,  one¬ 
dimensional  rules  can  be  thought  of  as  sequences.  Denoting  the  exact  integral  /,  (3.1) 
can  be  written 

The  theoretical  underpinnings  of  polynomial  methods  differ  substantially  from 
those  of  Riemann  sums.  Polynomials  are  fit  to  the  integrand,  and  the  resulting 
polynomial  function  is  then  integrated.  In  general,  the  nodes  are  the  roots  of  the 
fitted  polynomials.  For  this  reason,  these  methods  perform  exceptionally  well  on  very 
smooth  integrands  that  can  be  well  approximated  by  a  polynomial,  and  correspond¬ 
ingly  poorly  on  less  smooth,  badly  behaved  functions.  Three  different  polynomial 
methods  are  employed  here. 
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Gaussian  Quadrature 


Gaussian  quadrature  [6,15,17,36]  is  perhaps  the  best  known  polynomial-based  method, 
famed  for  its  2 n  —  1  degree  of  polynomial  exactness,  meaning  that  with  n  function 
evaluations,  all  integrands  which  are  polynomials  of  order  2 n  —  1  or  less  will  be  in¬ 
tegrated  exactly.  Gaussian  quadrature  utilizes  orthogonal  polynomials,  for  example 
Legendre,  Hermite,  or  Laguerre  polynomials.  Each  set  of  orthogonal  polynomials 
is  orthogonal  with  respect  to  a  different  inner  product,  the  kernel  of  which  is  the 
weighting  function  of  the  integral  approximated  by  the  associated  rule.  Legendre 
polynomials  are  orthogonal  with  respect  to  a  constant  weighting  function: 


J  Ll{x)L3{x)dx  = 


0  for  i  ^  j 

2i+I  for  ^  ~  j 


f{x)dx 


Also  of  importance  are  the  Hermite  polynomials,  which  are  orthogonal  with  respect 
to  the  Gaussian  weighting  function: 


_ x 2 /2 

~^Hl  (. x)Hj(x)dx 


0  for  i  ^  j 

n\- s/27t  for  i—j 


e-x2/2 

~7^!(x)dx 


S>.W)- 


When  approximating  expectations  of  functions  of  random  variables  (2.2),  the  choice 
of  orthogonal  polynomial  should  reflect  the  type  of  random  variable,  because  a  given 
probability  density  function  will  often  be  identical  to  some  orthogonal  polynomial’s 
weighting  function.  Obviously,  Legendre  polynomials  should  be  used  for  uniform 
random  variables,  and  Hermite  polynomials  for  normal  random  variables. 


Clenshaw-Curtis  Quadrature 

The  nodes  of  Clenshaw-Curtis  [4]  quadratures  are  at  the  roots  of  Chebyshev  poly¬ 
nomials.  In  the  past,  the  ease  with  which  they  could  be  computed  was  considered 
a  virtue  of  Clenshaw-Curtis  quadratures;  with  modern  computing  power,  this  is  no 
longer  a  significant  advantage.  Clenshaw-Curtis  only  achieves  n  polynomial  exact¬ 
ness,  but  is  known  to  exceed  its  theoretical  performance  expectations  [41].  A  further 
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advantage  is  nesting:  the  set  of  nodes  for  each  level  are  contained  in  the  nodes  for  all 
successive  levels.  In  section  3.1.2,  it  will  be  seen  that  this  is  desirable  for  the  sparse 
and  dimension-adaptive  grid  formulations,  where  nesting  allows  function  evaluations 
to  be  reused. 

Gauss-Kronrod-Patterson  Quadrature 

Gauss-Kronrod-Patterson  quadratures  [10,31,33]  are  an  extension  of  Gaussian  rules 
first  developed  by  Kronrod  and  then  iterated  by  Patterson  to  yield  a  sequence  of 
nested  quadratures.  The  ng  point  Gaussian  quadrature  base  level  has  a  polynomial 
exactness  of  2 ng  —  1,  and  each  successive  level  with  m  new  points  increases  the 
exactness  by  m.  The  overall  exactness  for  an  ng  point  Gaussian  quadrature  with 
X]  nr  Kronrod-Patterson  extension  points  is  then  2 ng  —  1  +  J2  m- 

In  section  4.3  we  find  that,  when  using  Gauss-Kronrod-Patterson  quadrature  be¬ 
ginning  with  the  three  point  Gaussian  rule,  shifting  the  index  by  one  yields  signifi¬ 
cantly  better  results  than  what  is  achieved  with  conventional  indexing,  that  is,  with 
the  first  level  corresponding  to  a  single  point.  We  find  a  similar  but  smaller  improve¬ 
ment  when  the  same  shift  is  applied  using  Clenshaw-Curtis  quadratures. 

3.1.2  In  Multiple  Dimensions 

One-dimensional  numerical  integration  schemes  are  extended  to  multiple  dimensions 
by  tensor  products,  which,  in  all  implementations  considered  here,  are  more  plainly 
Cartesian  products.  For  example,  suppose  <3>io  is  the  ten  point  Legendre  polynomial 
Gaussian  quadrature.  Then  <hi0  ®  $10  is  the  one  hundred  point  two-dimensional  inte¬ 
gral  approximation,  and  $10  $10  $10  is  the  one  thousand  point  three-dimensional 

integral  approximation.  It  is  through  this  exponential  increase  in  points  necessary  to 
maintain  a  constant  level  of  discretization  with  increasing  dimension  that  the  Curse 
of  Dimensionality  enters  numerical  integration.  This  example  is  a  ’full  grid’  multidi¬ 
mensional  scheme.  In  the  coming  sections,  sparse  grids,  dimension-adaptivity,  and, 
the  alternatives  to  Cartesian  product  methods  in  high  dimensions,  Monte  Carlo  and 
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quasi-Monte  Carlo  will  be  discussed. 


Full  Grids 

The  d-dimensional,  level  q  full  grid  collocation  integral  approximation,  which  we  will 
denote  Fq,  is  the  tensor  product  of  d  one  dimensional,  level  q  quadratures: 


Fqd[f]  =  <!>q®---®%[f} 

=  Eqli  •  •  •  EjLi  wilt...,idf(x  n, xid)  (3-2) 

Wilt...,id  =wh-  ...■  wid, 

where  rnq  is  the  number  of  points  in  Qq.  The  full  grid  formulation  is  written  here  as 
having  identical  component  quadratures  in  each  dimension  for  concision,  but  this  is 
not  a  requirement:  any  level  of  any  valid  rule  may  be  used.  This  is  often  advantageous 
to  keep  in  mind;  a  function  that  is  poorly  behaved  in  some  dimensions  and  relatively 
smooth  in  others  can  be  more  efficiently  integrated  by  using  high  level  approximations 
in  the  troublesome  dimensions  and  lower  level  approximations  in  the  others,  rather 
than  by  using  the  same  high  level  quadrature  in  every  dimension.  The  dimension- 
adaptive  introduced  in  section  3.1.2  does  this  automatically. 

With  some  information  about  the  function  being  integrated,  the  error  for  full 
grid  collocation  with  polynomial  based  component  quadratures  can  be  bounded.  For 
functions  in  Cr  (functions  with  bounded  mixed  derivatives  up  to  order  r),  the  error 
of  full  grid  collocation  is  estimated  by 

m  -  F^f}\  =  (3.3) 

where  nq  is  the  total  number  of  points  in  the  grid.  The  factor  of  d  in  the  exponent 
is  a  manifestation  of  ’The  Curse  of  Dimensionality’.  As  dimension  grows,  it  causes 
the  error  order  to  grow  as  well  at  an  exponential  rate.  For  low  dimensional,  smooth 
integrands,  a  full  grid  is  often  the  best  choice  among  non-adaptive  methods. 
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Sparse  Grids 


In  1963,  Smolyak  introduced  what  is  now  known  as  Smolyak’s  formula  and  is  the 
underlying  formulation  of  all  sparse  grid  methods  [12,21,28,34,38,43].  Rather  than 
using  one  high  order  tensor  product,  the  sparse  grid  is  a  sum  of  lower  order  Cartesian 
products,  and  tends  to  perform  better  in  slightly  higher  dimension  than  full  grid 
approaches.  Define  the  level  %  difference  quadrature 

A;  =  A!  =  $!. 

Smolyak’s  formula  can  then  be  written 

5,1/]=  V  (3.4) 

\i\i<q+d 

where  =  A  +  •  •  •  +  id- 

At  this  point  it  becomes  clear  why  nested  component  quadratures  are  preferred 
for  sparse  grids.  Nesting  causes  points  of  different  grids  in  the  sum  to  coincide,  and 
the  number  of  common  points  increases  with  both  the  level  and  dimension  of  the 
sparse  grid. 

The  error  for  a  sparse  grid  with  polynomial  based  component  quadratures  can  be 
estimated  in  the  same  fashion  as  the  full  grid  case: 

W}  -  S,1/]|  =  0(n~rlog(nq)^+»).  (3.5) 

As  with  full  grid  collocation,  the  factor  d  eventually  causes  sparse  grids  to  fall  to  ’The 
Curse’  as  well,  but  not  as  quickly.  The  penalty  incurred  by  dimension  is  less  severe 
because  of  the  log(nq)  term,  and  consequently  sparse  grid  collocation  is  often  superior 
in  mid- dimensional  situations.  Fig.  3-1  shows  full  and  sparse  grid  examples. 


Full  grid  Sparse  grid 


Figure  3-1:  Left:  Full  grid  with  Legendre  polynomial  Gaussian  quadrature.  Right: 
Sparse  grid  with  Gauss-Kronrod-Patterson  quadrature. 

Dimension- Adaptivity 

Reference  [13]  presents  an  adaptive  algorithm  which  places  points  according  to  online 
estimation  of  its  convergence.  The  basic  building  block  of  the  dimension-adaptive 
algorithm  is  the  difference  grid  A  defined  in  the  previous  section.  Let  E.,  [/]  =  /[/]  — 
Q J/]  be  the  error  of  some  quadrature  Qj,  be  it  one  or  multidimensional.  Because  the 
error  is  converging  to  zero,  E,  >>  E,+1,  and  E,  ~  E,  —  Ei+1.  Then 

Ej  «  Ej  —  Ei+1  =  I  —  Qj  —  I  +  Qj+i  =  Qi+i  —  Q,  =  Ai+1 . 

The  algorithm  builds  an  approximation  from  difference  grids,  and  their  magnitude, 
which  estimates  the  convergence  of  a  dimension  or  combination  of  dimensions,  indi¬ 
cates  which  difference  grids  should  be  evaluated  next.  For  example,  suppose  that  for 
a  two  dimensional  integrand  the  evaluated  difference  grid  with  the  largest  magnitude 
is  A28  A3.  Then,  if  they  are  valid  difference  grids,  the  next  two  grids  to  be  evaluated 
will  be  A3  <8)  A3  and  A2  8  A4.  Fig.  3-2  illustrates  the  behavior  of  the  algorithm 
on  a  test  integral  from  section  4.1.  To  demonstrate  the  versatility  of  the  algorithm, 
the  space  has  been  divided  into  nine  identical  elements,  inside  each  of  which  the 
dimension-adaptive  algorithm  was  executed.  Multi-elements  [20]  is  another  powerful 
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approach  to  numerical  integration,  but  which  is  not  explored  here.  The  function  is 
C°  along  the  central  axes  of  the  space,  and  analytic  elsewhere.  Grids  that  are  of  high 
order  in  the  dimension  perpendicular  to  the  C°  ridges  are  evaluated,  while  relatively 
low  order  grids  are  computed  where  the  function  is  better  behaved. 


Evaluations 


Figure  3-2:  Nine  element  Clenshaw-Curtis  dimension-adaptive  integration  of  the  two 
dimensional  instantiation  of  test  function  five  from  [11],  The  nodes  chosen  by  the 
algorithm  are  shown  on  the  left,  and  the  function  evaluations  at  those  points  on  the 
right. 


The  algorithm  (from  [13])  is  as  follows: 
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Algorithm: 

DimAdapt(f) 

i=[l,...,l] 

0  =  0 
A  =  {i} 

r  =  A  j/ 

V  =  9i 

while  rj  >  to/ 

im  =  {i  G  A  :  cjfj  =  max  g} 

A  =  A\im 
O  =  O  U  im 

1  =  11-9^ 
for  k  —  1  :  d 

j  bn.  "h 

if  j  —  eq  G  O  for  all  q  —  1  :  d 
A  =  A  U  j 

r  =  r  +  Aj/ 

V  =  V  +  flj 
end  if 
end  for 
end  while 
return  r 


Symbols: 

O  old  index  set 
A  active  index  set 
A •/  integral  increment 

g*  local  error  indicator 
global  error  estimate 
kth  unit  vector 
error  tolerance 
computed  integral  value 


n 

V 

tol 

r 


Monte  Carlo  and  Quasi-Monte  Carlo 


For  all  but  trivial  functions,  Cartesian  product  based  methods  will  invariably  be¬ 
come  ineffective  in  high  enough  dimension.  Although  it  is  not  the  focus  of  this  work, 
it  is  worth  discussing  the  high- dimensional  alternative:  Monte  Carlo  methods.  Pure 
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Monte  Carlo,  or  pseudo-randomness,  converges  according  to  the  central  limit  theorem 
at  1/y/n,  independently  of  dimension.  Monte  Carlo  methods  also  do  not  depend  on 
the  smoothness  of  the  integrand,  and  for  badly  behaved  functions  are  sometimes  supe¬ 
rior  even  in  low  dimensions.  Quasi- Monte  Carlo  [3,27,39]  methods  offer  improvement 
to  1/n  optimally  and  log (n)d/n  in  the  worst  case.  None  of  these  are  attractive  rates, 
but  in  high  dimensions  it  is  often  the  best  there  is.  Fig.  3-3  shows  pseudo-random 
and  Sobol  sequence  quasi-random  points. 


Pseudo-random 


Figure  3-3:  Left:  100  pseudo-random  points.  Right:  100  quasi-random  points  from 
the  Sobol  sequence. 


3.2  Interpolation 

Interpolation  [2,7,21]  is  closely  tied  to  numerical  integration,  and  often  utilizes  the 
same  function  samples.  In  fact,  with  the  exception  of  Monte  Carlo,  all  of  the  al¬ 
gorithms  discussed  in  this  chapter  can  be  used  to  interpolate  as  well  as  integrate 
functions  (Monte  Carlo  is  in  a  sense  a  more  direct  but  sometimes  less  efficient  route 
to  what  interpolation  provides).  Essentially,  rather  than  fitting  polynomials  to  a 
function  and  integrating  the  fit,  we  are  now  evaluating  it. 

The  one- dimensional  barycentric  interpolation  formula  of  the  second  form  is  prefer- 
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able  for  its  numerical  stability  and  computationally  efficiency.  It  is  given  by 


wiw  = 


m  f(„\ 
1  x—Xi’'  1  D 

\r^mq  Wj 

£—‘i=  1  x—Xi 


Wi  =  n — 

j=i  ^ i 

j¥=i 


(3.6) 


Full,  sparse,  and  dimension-adaptive  interpola.nts  formulas  are  constructed  by  tak¬ 
ing  Cartesian  products  of  one  dimensional  interpolants,  as  with  numerical  integration. 

The  intent  is  to  use  interpolation  to  generate  PDF’s,  i.e.  histograms.  A  direct  way 
to  accomplish  this  without  interpolation  is  with  Monte  Carlo:  evaluate  a  function  or 
simulate  a  system  at  a  bunch  of  random  points  and  sort  the  outcomes  into  bins.  So 
why  then  is  interpolation  necessary?  Rather  than  direct  evaluation,  which  can  be 
costly,  instead  apply  Monte  Carlo  to  the  interpolant,  which  is  always  cheap  to  eval¬ 
uate,  and  generate  histograms  of  interpolated  function  values  or  system  simulations. 
Again  consider  the  electric  ship;  one  evaluation  takes  ten  minutes,  so  directly  applying 
Monte  Carlo  to  create  a  1,000  point  histogram  would  be  time  consuming.  However, 
suppose  the  response  can  be  characterized  with  fair  accuracy  using  100  collocation 
points.  The  simulation  can  now  be  interpolated  1,000  times  to  produce  nearly  the 
same  histogram  (the  similarity  increases  with  the  number  of  collocation  points  used) 
in  a  trivial  length  of  time. 

Again  consider  the  double  pendulum  example  of  section  2.2.  Fig.  3-4  shows 
interpolated  trajectories  overlaid  on  trajectories  obtained  by  direct  simulation  corre¬ 
sponding  to  the  same  realization  of  the  uncertain  parameters.  Fig.  3-5  shows  two 
sets  of  histograms  generated  by  the  same  Monte  Carlo  points,  on  the  left  generated 
by  interpolation  and  on  the  right  by  direct  simulation. 

There  is  strong  agreement  between  the  simulated  and  interpolated  trajectories;  it 
is  only  near  the  end  when  the  variance  has  grown  large  that  any  disparity  is  noticeable. 

The  histograms  also  agree  fairly  well.  All  have  multiple  modes,  meaning  that 
each  arm  can  have  substantially  different  positions  and  velocities  depending  on  where 
in  the  random  space  the  system  is  realized,  and  the  interpolated  histograms  have 
captured  this  feature. 
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15 


Directly  evaluated  and  interpolated  trajectories 


Actual 


Figure  3-4:  Interpolated  and  directly  evaluated  trajectories  of  double  pendulum 
states. 


3.3  Summary 

An  overview  of  collocation  algorithms  for  numerically  solve  stochastic  differential 
equations  is  given,  along  procedures  for  extracting  statistical  moments  and  PDF’s 
from  the  resulting  solutions.  Guidelines  for  using  each  algorithm  are  suggested,  specif¬ 
ically  that  full  and  sparse  grids  perform  well  on  smooth  functions  in  low  and  mid 
dimensions,  and  that  Monte  Carlo  or  quasi-Monte  Carlo  will  eventually  become  su- 
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Interpolated  histograms 


-2  0  2 


Directly  evaluated  histograms 


-2  0  2 


Figure  3-5:  Interpolated  and  directly  evaluated  histograms  of  double  pendulum  states 
at  2.5  seconds. 

perior  with  increasing  dimension,  or  even  in  lower  dimensions  if  the  function  is  poorly 
behaved.  The  dimension- adaptive  algorithm  from  [13]  in  many  cases  will  outperform 
the  non-adaptive  collocation  methods. 
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Chapter  4 


Comparison  of  Integration 
Methods 


Numerical  integration  is  perhaps  the  simplest  setting  in  which  to  compare  the  per¬ 
formance  of  the  algorithms  of  the  previous  chapter.  In  this  chapter,  a  comparison  is 
given  using  the  scalar  test  functions  from  [11],  In  section  4.3,  numerical  results  and 
theoretical  justification  are  presented  for  a  modification  to  conventional  sparse  grids 
yielding  improved  performance. 


4.1  Test  functions 


The  six  test  functions  from  [11]  are: 


OSCILLATORY  /(x) 
PRODUCT  PEAK  /(x) 
CORNER  PEAK  /(x) 
GAUSSIAN  /(x) 

C°  /(x) 

DISCONTINUOUS  /(x) 


COs(27TUi  +  J2i= 1  UiXi) 

nf=iKr2  +  (xi  -Ui)2)-1 

(1  +  Eti  a^)-(d+1) 

exp(—  Eli  a^Xi-Ui)2) 
exp(—  Ei=i  cn\xi  -  Ui\ ) 

{0  x\  >  u\  or  X2  >  ii2 

exp(Eti  aixi)  otherwise 
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where  a  and  u  are  tunable  parameters  which  respectively  determine  the  difficulty  of 
the  functions  to  integrate  and  the  location  of  their  features  in  the  region  of  integration. 
Plots  of  the  two  dimensional  instantiation  of  each  function  are  shown  in  Fig.  4.1.  In 


Oscillatory 


Product  Peak 


C° 


Gaussian 


Figure  4-1:  Two  dimensional  test  functions 


the  numerical  experiments  of  the  next  section,  a  was  chosen  to  make  each  integral 
roughly  the  same  difficulty,  and  u  was  generated  randomly  over  a  set  of  trials. 
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4.2  Convergence  of  Methods 


The  theoretical  error  estimates  from  section  3.1.2  say  that  full  grids  will  perform 
best  on  smooth  integrands  in  low  dimensions,  sparse  grids  on  smooth  integrals  in 
mid- dimensional  ranges,  and  Monte  Carlo  or  quasi-Monte  Carlo  on  high  dimensional 
integrands  and  even  non-smooth  integrands  in  low  dimensions.  Note  that  the  compar¬ 
ison  given  here  is  summary,  but  adequate;  it  is  not  necessary  to  show  the  performance 
of  every  combination  of  one  dimensional  quadrature  and  multidimensional  algorithm 
on  each  test  function  in  a  full  range  of  dimensions  to  capture  the  essence  of  when 
each  method  is  appropriate. 

4.2.1  Non-adaptive  methods 

Fig.  4.2.1  shows  convergence  plots  comparing  full  grid  with  Legendre  polynomial 
Gaussian  Quadrature,  sparse  grid  with  Clenshaw-Curtis,  Monte  Carlo,  and  quasi- 
Monte  Carlo  on  the  oscillatory,  product  peak,  and  corner  peak  test  functions.  Fig. 
4.2.1  shows  the  same  for  the  Gaussian,  C°,  and  discontinuous  test  functions. 

The  oscillatory  function  is  analytic,  and  as  would  be  expected  the  collocation 
algorithms  are  superior  to  Monte  Carlo  and  quasi-Monte  Carlo  even  into  ten  dimen¬ 
sions,  although  the  performance  gap  is  decreasing.  The  full  grid  is  superior  in  two 
dimensions,  roughly  equivalent  to  the  sparse  grid  in  six  dimensions,  and  inferior  in 
ten,  also  what  is  expected  from  the  theoretical  bounds. 

The  product  peak,  corner  peak,  and  Gaussian  functions  have  similar  stories:  full 
grids  are  best  in  low  dimensions,  sparse  grids  in  mid  dimensions,  and  Monte  Carlo 
and  quasi-Monte  Carlo,  if  not  the  favorite  by  dimension  ten,  soon  will  be.  Notice  the 
corner  peak  plot  for  ten  dimensions:  Monte  Carlo  is  dramatically  more  efficient  than 
quasi-Monte  Carlo.  Quasi- Monte  Carlo  has  1/iV  convergence  in  the  optimal  case,  but 
depends  on  dimension  in  the  worst  case,  as  is  clearly  happening  here. 

The  continuous  and  discontinuous  functions  are  considerably  more  difficult  than 
the  rest.  Both  have  very  low  smoothnesses,  and  are  therefore  difficult  to  fit  poly¬ 
nomials  to.  Monte  Carlo  and  quasi-Monte  Carlo  are  arguably  for  the  continuous 
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2  dimensions 


6  dimensions 


10  dimensions 


10 


Figure  4-2:  Log-log  plots  of  error  versus  number  of  function  evaluations  for  product 
peak,  oscillatory,  and  corner  peak  test  functions.  Legend  entries  correspond  to:  FG- 
GQ  -  Full  grid  with  Gaussian  Quadrature,  SG-CC  -  Sparse  grid  with  Clenshaw- Curtis, 
MC  -  Monte  Carlo,  QMC  -  Quasi-Monte  Carlo. 
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Figure  4-3:  Log-log  plots  of  error  versus  number  of  function  evaluations  for  Gaussian, 
(7°,  and  discontinuous  test  functions. 


function  and  clearly  for  the  discontinuous  function  the  best  choice  even  as  low  as  two 
dimensions. 
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4.2.2  Dimension-adaptivity 

Figures  4.2.2  and  4.2.2  compare  full  and  sparse  grids  with  the  dimension-adaptive 
algorithm  discussed  in  section  3.1.2,  [13].  The  dimension- adaptive  algorithm  is  as 
good  as  or  superior  to  the  best  of  the  other  methods  in  nearly  all  cases,  with  the 
adaptivity  ’running  astray’  only  on  the  discontinuous  function,  for  which  a  Monte 
Carlo  would  be  superior. 


Number  of  function  evaluations,  log  scale 


Figure  4-4:  Log-log  plots  of  error  versus  number  of  function  evaluations.  Sparse  grid 
with  Gauss-Kronrod-Patterson  (SG-KP),  full  grid  with  Gaussian  Quadrature  (FG- 
GQ),  and  dimension-adaptivity  with  Gauss-Kronrod-Patterson  (DA-KP)  in  three  di¬ 
mensions  are  shown. 


The  performance  of  each  method  varies  significantly  across  the  different  test  func¬ 
tions.  No  one  method  is  universally  superior,  but  prior  knowledge  of  the  integrand, 
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Oscillatory 


Product  Peak 


Continuous  Discontinuous 


Number  of  function  evaluations,  log  scale 


Figure  4-5:  Log-log  plots  of  error  versus  number  of  function  evaluations  for  algorithms 
in  six  dimensions  are  shown. 


such  as  smoothness  and  dimension,  can  provide  information  that  is  highly  relevant  to 
choosing  an  approach.  It  is  observed  that  among  the  non-adaptive  methods,  full  grids 
are  superior  in  low  dimensions,  sparse  grids  in  mid  dimensions,  and  Monte  Carlo  and 
quasi-Monte  Carlo  in  high  dimensions  or  for  non-smooth  integrands.  The  dimension- 
adaptive  algorithm  in  most  cases  outperforms  the  other  grid  methods,  but  should  be 
used  cautiously,  particularly  on  badly  behaved  functions. 
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4.3  Sparse  Grids  with  Shifted  Indexing 


A  slight  modification  to  the  indexing  of  a  sparse  grid  yields  significant  gains  in  conver¬ 
gence  in  mid  dimensional  ranges,  particularly  for  smoother  functions.  Over  the  next 
few  sections,  the  modification  is  described,  then  theoretically  and  computationally 
bounded,  and  lastly  results  from  numerical  experiments  are  shown. 


4.3.1  One-Dimensional  Quadratures  and  Shifting 

When  using  Kronrod-Patterson  quadrature  beginning  with  the  three  point  Gaussian 
rule,  shifting  the  index  by  one  yields  significantly  better  results  than  what  is  achieved 
with  conventional  indexing,  that  is,  with  the  first  level  corresponding  to  only  one 
point.  A  similar  but  smaller  improvement  is  found  when  the  same  shift  is  applied 
using  Clenshaw- Curtis  quadratures.  A  sparse  grid  with  indices  shifted  a  times  will 
be  denoted  S^a  and  the  corresponding  one- dimensional  basis  sequence  by  <f>". 


4.3.2  Error  Bounds 

Effect  of  the  Shift  on  Classical  Bounds 

The  error  of  the  shifted  grid  can  be  bounded  in  the  same  fashion  as  a  conventionally 
indexed  sparse  grid,  by  allowing  the  shift  to  propagate  through  to  the  final  estimate. 
Following  the  proof  of  [43]  (see  also  [28]),  the  derivation  of  an  error  bound  for  the 
d-dimensional,  once-shifted  sparse  grid  Sf{ ,  is  sketched.  Begin  by  estimating  the  error 
of  a  one-dimensional  quadrature  as 


||/i  -$i||  <  cr  ■2“ri,  i  >  0  (4.1) 

where  ||-||  =  max  |  •  |,  <f>0[/]  =  0,  ||/i||  =  2,  and  /  has  bounded  mixed  derivatives  up  to 
order  r.  The  result  of  removing  the  lowest  one-dimensional  quadrature  is  a  shifting 
of  indices;  the  new  error  estimate  is 


h  -  $  •  <  cr  ■  2~r(i+1)  i>  0. 


(4.2) 
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The  magnitudes  of  the  one- dimensional  once-shifted  difference  grids  are  bounded  by 


K  <  h~K  +  <cr.2-’'<i‘+1».(l  +  2-).  (4.3) 

Making  use  of  the  property  ||kl  ®  £>||  =  ||kL||  ■  ||5||,  we  have 


- s,% l<  E  IK  ■■■■■  K  ■ 

l*ll<9 

h  -  tj+i_|4  +  2  •  Id-  s£,  ,  so  that  (4.4) 


IWl-StbOf/-1^-^),  (4.5) 

which  is  very  similar  to  Smolyak’s  original  estimate: 

\ld[f]  ~  -S', 1 1  =  0(qd~l  ■  2~rq).  (4.6) 

The  new  bound  differs  by  a  factor  of  2~rd  from  the  original,  which  suggests  an  advan¬ 
tage  increasing  with  dimension.  If  the  one-dimensional  quadrature  indexing  is  shifted 
by  a  rather  than  1,  a  substitution  in  the  above  bound  gives 

h[f]-Sdq,«  ~  0(qd~1  ■  2~r(q+ad)). 
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Next,  a  point-based  bound  is  obtained  for  the  once-shifted  case,  by  estimating 
the  total  number  of  points  in  the  sparse  grid  and  substituting  for  q.  Assuming  the 
number  of  points  in  each  one-dimensional  quadrature  rrit  <  2\  the  total  number  of 
points  n  is  approximated  by 

n  <  2q+d  •  j  j  =  0(qd~1  ■  2q+d). 

Substituting  this  into  the  level  bound  above  gives 


The  inequality  2q+d  <  c^n,  q  <  log (c^n)  —  d  is  then  substituted  to  arrive  at: 

\ld[f]  -  ^i||  =  0{n-r  ■  (log(n)  -  d)^1^). 

The  point-based  bound  corresponding  to  (4.6),  the  non-shifted  estimate,  is 

\ld[f]  -  Sd\\  =  0{n~r  ■  (log(n))(d-1)(r+1)), 

which  again  differs  by  a  term  containing  d,  although  here  it  is  in  the  base  instead  of 
the  exponent. 


The  estimates  in  (4.1)  and  (4.2)  are  very  crude,  and  as  a  result  the  bound  in 
(4.5)  fails  to  thoroughly  capture  the  behavior  of  the  shifted  scheme  seen  in  numeri¬ 
cal  experiments,  namely  that  there  is  an  improvement  in  intermediate  but  not  high 
dimensions.  In  the  next  two  sections,  an  exact  calculation  for  the  number  of  points 
and  a  more  precise  formulation  of  one-dimensional  quadrature  errors  is  given,  and 
used  them  with  (4.3)  and  (4.4)  to  arrive  at  a  more  accurate  description  of  the  shifted 
sparse  grid’s  performance. 
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The  Number  of  Points 


The  exact  number  of  points  in  the  nested  sparse  grid  S^a  can  be  directly  computed 
by  Theorem  1  in  [34],  Define  Si  as  the  number  of  points  used  in  <f>f  and  not  4>f_1  and 
set  k  =  q  —  d.  The  total  number  of  points  is  then  given  by 


n(d  +  k,d)  =  chk 

3=0  W 


where  (^j  =  0  if  j  >  d,  and  c  is  defined  by 


cj,k  5Z 


1-3  + 1  X 

^  -r-  1  ,  1  —  1  h 

j  c  —  l,fc— VI  J  -L,  .  .  .  ,  A/, 

1  ^0 


Cq  ,fc  1- 


One-Dimensional  Bounds  through  Peano’s  Theorem 

The  Peano  kernel  [6]  of  a  quadrature  is  defined  as 


Kr(t)  =  j,  ■  E,[(x 

-tn 

[  x  —  t 

X  >  t 

x-t)+  =  < 

l 0 

X  <  t 

where  p  is  the  polynomial  exactness  and  Et  the  error  of  the  quadrature  <f>j  (E,  [g\  =  0 
for  all  polynomials  g  of  order  less  than  or  equal  to  p).  The  quadrature  error  can  be 
written  as 

E,\f]  =  /[/]  -  $,[/]  =  j Kp(t)f<r+'\t)dt. 

for  all  /  in  Cp+1.  The  error  can  then  be  approximated  by  either  of 

Ei[f]  <  \\KP\\  j \f^)(x)\dx  or  Ei[f]  <  ||/(p+1)(t)||  j  \Kp\dx 
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where  ||-||  denotes  the  sup  norm.  The  magnitude  of  the  difference  grid  (4.3)  is  esti¬ 
mated  with  the  second  inequality  above  by  fixing  /(p+1')(a;)  at  one  for  all  p.  This 
constraint  on  the  function  derivatives  is  a  serious  simplification,  but  can  be  relaxed 
for  cases  where  some  properties  of  the  function  are  known. 


4.3.3  Numerical  Results 

As  shown  below,  in  the  non-shifted  case  Kronrod-Patterson  generally  performs  marginally 
better  than  Clenshaw-Curtis  based  sparse  grids.  Kronrod-Patterson’s  polynomial 
exactness  is  only  slightly  superior  to  that  of  Clenshaw-Curtis,  and  furthermore, 
Clenshaw-Curtis  has  been  known  to  exceed  its  theoretical  expectations  [41].  Dif¬ 
ferences  in  polynomial  exactness  are  mostly  due  to  the  three-point  levels  of  each 
sequence:  the  three  point  Kronrod-Patterson  rule  has  Gaussian  quadrature’s  poly¬ 
nomial  exactness  of  five,  while  Clenshaw  Curtis  only  has  a  polynomial  exactness  of 
two. 

Computed  Error  Bounds 

Now,  using  (4.4)  and  the  results  of  the  sections  4.3.2  and  4.3.2,  bounds  for  the  con¬ 
ventional  and  shifted  grids  can  be  realistically  computed  and  compared  (Fig.  4-6). 
Firstly,  observe  that  the  rate  of  convergence  of  both  shifted  grids  is  superior  in  all 
dimensions,  and,  according  to  the  crude  estimate  (4.5),  further  distancing  itself  from 
the  non-shifted  grid  as  the  dimensions  grows.  In  the  Kronrod-Patterson  plots  for  a 
single  shift,  a  ’’knee”  occurs  in  the  non-shifted  curves,  after  which  a  similar  slope  is 
achieved.  Assessment  of  performance  for  the  two  rules  with  and  without  the  shift  is 
best  summarized  by  noting  where  the  shifted  and  non-shifted  curves  cross,  i.e.  where 
the  accuracy  per  effort  is  identical.  By  this  criteria,  it  is  evident  that  in  higher  di¬ 
mensions  a  desired  accuracy  may  be  attainable  with  a  non-shifted  grid  before  the 
shifted  grid  has  become  effective.  The  Kronrod-Patterson  case,  however,  shows  that 
the  benefits  of  the  shift  can  be  realized  at  much  more  relaxed  tolerances. 

Another  point  worth  considering  is  the  performance  of  sparse  grids  shifted  more 
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Three  dimension 


Five  dimensions 


Seven  dimensions 


Nine  dimensions 


Non-shifted 

Once-shifted 


O  Twice-shifted 


Figure  4-6:  Theoretical  error  convergence  of  conventional,  once-shifted,  and  twice- 
shifted  sparse  grids  with  Clenshaw-Curtis  and  Kronrod-Patterson  basis  quadratures. 

than  once.  Included  in  Fig.  4-6  are  theoretical  error  curves  for  twice-shifted  sparse 
grids.  The  rates  of  convergence  are  even  higher,  and  for  Clenshaw-Curtis  the  crossing 
with  the  non-shifted  curve  occurs  at  approximately  the  same  location  as  that  of  the 
once-shifted  curve.  However,  the  dramatic  increase  in  points  associated  with  each 
shift  can  make  higher  order  shifted  grids  unwieldly  even  in  moderate  dimensions,  and 
we  do  not  have  numerical  experiments  to  verify  these  predictions. 

Numerical  Examples 

As  in  section  4.2,  the  performance  of  sparse  grids  with  shifted  indexing  is  evaluated 
on  test  integrals  from  [11],  specifically  the  oscillatory,  product  peak,  corner  peak,  and 
Gaussian  functions  (the  smooth  ones). 

Fig.  4-7  shows  the  performance  of  regular  and  shifted  sparse  grids.  For  Kronrod- 
Patterson  based  grids,  significant  improvement  is  seen  for  the  shown  test  integrals, 
with  the  performance  improvement  due  to  the  shift  most  noteworthy  at  dimension 
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four.  As  predicted,  only  a  minor  improvement  is  evident  with  the  shift  only  in  low 
dimensions  when  using  Clenshaw- Curtis;  at  high  dimensions,  shifting  is  clearly  un¬ 
suitable. 


Two  dimensions  Four  dimensions  Six  dimensions 


101  102  103  101  102  103  104  101  102  103  104  105 


Figure  4-7:  Error  versus  function  evaluations  for  regular  and  shifted  sparse  grids  using 
Clenshaw- Curtis  and  Kronrod-Patterson  quadratures 


4.3.4  Discussion 

It  has  been  shown  that  for  smooth  functions,  shifting  the  indexing  of  a  sparse  grid 
yields  improved  accuracy  which  sometimes  outweighs  the  associated  increase  in  points. 
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The  modified  scheme  is  more  effective  with  Kronrod-Patterson  than  with  Clenshaw- 
Curtis  sequences,  a  difference  which  is  suggested  by  the  polynomial  exactnesses  of  the 
three  point  levels  of  each  sequence.  Intuitively,  shifting  succeeds  because  high  fidelity 
on  one  hyperplane  is  traded  for  evenly  distributed  high  fidelity  on  three  hyperplanes. 
Hence,  one  would  expect  inferior  performance  for  functions  with  weak  dimensional 
coupling. 

4.4  Summary 

The  algorithms  of  the  previous  chapter  are  compared  on  six  standard  test  functions, 
and  expectations  for  the  performance  of  each  algorithm  were  confirmed.  A  mod¬ 
ification  to  the  conventional  sparse  grid  algorithm  was  presented  and  justified  via 
theoretical  bounds  and  numerical  examples. 
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Chapter  5 


Stochastic  Simulation  of  the 
All-Electric  Ship 


In  this  chapter,  collocation  algorithms  are  applied  to  three  large-scale  models,  two 
describing  an  electric  ship  and  the  other  a  pulse  power  system.  Sensitivity  is  inferred 
for  each  system  to  demonstrate  how  information  from  stochastic  simulation  can  be 
interpreted. 


5.1  ONR  IPS  Testbed 

5.1.1  Deterministic  Simulation 

We  simulate  the  Office  of  Naval  Research  Integrated  Power  System  testbed  (ONR 
IPS)  [23,32]  (Fig.  5-1),  varying  parameters  from  the  port  and  starboard  AC  systems, 
which  contains  the  starboard  generator,  bus,  and  propulsion. 

First,  we  examine  what  happens  in  the  model  through  deterministic  simulation. 
In  zone  three,  a  constant  power  load,  a  proportional-integral  controller,  and  two 
ship  service  converter  modules  are  turned  on,  all  at  time  zero.  The  zone  three  port 
voltage,  our  representative  voltage  for  zone  three,  rises  from  the  start  until  about 
0.15  seconds  when  it  reaches  400  volts.  At  time  0.1  seconds,  the  starboard  propulsion 
current  becomes  non  zero  -  this  is  when  the  motor  and  motor  control  turn  on.  At 
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PS:  Power  Supply  SSIM:  Ship  Service  Inverter  Module 

(3-<|)  AC  to  DC)  (DC-AC  Converter) 

SSCM:  Ship  Service  Conerter  Module  LB:  Local  Bank 

(DC-DC  Converter)  MC:  Motor  Controller 

CPL:  Constant  Power  Load 

Figure  5-1:  Schematic  for  ONR  IPS  Testbed. 


0.4  seconds,  the  system  is  loaded  from  the  ship  service  inverter  module  in  zone  one, 
and  the  power  supply  from  the  starboard  AC  system  increases  correspondingly  (Fig. 
5-2). 

5.1.2  Probabilistic  Simulation 

In  addition  to  determining  mean  behaviors,  multidimensional  techniques  can  be  em¬ 
ployed  for  sensitivity  analysis,  allowing  for  identification  of  not  only  individual  sen¬ 
sitive  parameters,  but  also  sensitive  combinations  of  parameters.  In  onr  simulations, 
we  use  both  the  integral  (IS)  and  maximum  (MS)  over  time  of  the  variance  solution 
as  a  sensitivity  metric. 

IS  —  J  var(x(t))dt,  MS  =  m ax{var(x(t))}  (5.1) 


54 


Zone  3  Port  Voltage 


Starboard  AC  System  Starboard  Power  Supply  Current  a  coordinate 


Figure  5-2:  Deterministic  trajectories  for  ONR  IPS. 

The  analysis  in  the  next  section  was  done  using  the  integral  metric. 

Non- Adaptive  Simulation 

We  analyze  results  from  a  three  dimensional,  level  four  full  grid  collocation  simulation 
varying  each  uncertain  parameter  ±1%  from  its  mean.  We  considered  all  triples  drawn 
from  the  set: 

•  Bus 

1,  7)  C  -  Shunt  Capacitance 
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2,  8)  Cf  -  Filter  Capacitance 

3,  9)  Lf  -  Filter  Inductance 

•  Propulsion 

4,  10)  C  -  Filter  Capacitance 

5,  11)  L  -  Filter  Inductance 

6,  12)  Induction  Motor  Mechanical  Load 

where  parameters  1-6  are  starboard  and  7-12  port.  Fig.  5-3  shows  time  solutions  for 
the  triple  {1,  2,  3}. 


Mean  solution  through  time  with  standard  deviation  error  bars:  0  to  0.5  seconds 


Seconds 


Figure  5-3:  Starboard  AC  system  generator  ahe-reference  frame  c  current  mean 
(curve)  and  standard  deviation  (error  bars)  solutions  for  the  triple  containing  star¬ 
board  bus  shunt  capacitance,  filter  capacitance  and  filter  inductance  for  ONR  IPS. 

There  are  many  ways  to  structure  results  from  this  sort  of  simulation,  for  example 
the  most  sensitive  triple  for  a  specific  output:  we  find  that  triple  {1,  2,  3}  is  in  fact  the 
most  sensitive  triple  for  the  starboard  generator  currents  according  to  the  integrated 
variance  metric,  and  {1,2,5}  by  the  maximum  metric.  We  present  results  that  lend 
themselves  easily  to  graphical  analysis,  but  they  are  not  the  only  perspectives. 
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For  simulations  in  n  dimensions,  the  interaction  of  pairs  of  uncertain  parameters 
can  be  viewed  by  considering  all  n-tuples  containing  the  same  n  —  2  parameters.  Fig. 
5-4  shows  all  triples  containing  the  starboard  bus  filter  inductance  (3).  It  can  be  seen 
that  the  most  influential  pair  on  this  plot  is  {1,2},  which  corresponds  to  what  we 
previously  identified  as  the  most  sensitive  triple. 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 
11 


Low  sensitivity  High  Sensitivity 


Relative  Integrated  Variance 


1  2  3  4  5  6  7  8  9  10  11  12 


Figure  5-4:  Two  dimensional  integrated  variance  analysis  of  the  starboard  AC  system 
generator  afec-reference  frame  c  current:  all  triples  containing  the  starboard  bus  filter 
inductance  (3)  for  ONR  IPS. 

Taking  the  most  influential  triple,  we  now  compare  the  sensitivity  of  different 
output  states.  Since  we  are  now  comparing  different  outputs  with  differing  units,  we 
normalize  the  integrated  variance  by  the  integral  of  the  square  of  the  deterministic 
solution.  The  output  states  corresponding  to  each  number  are  listed  below. 

•  Zone  3  DC  states 

1)  Load  Voltage 

2)  Starboard  Voltage 

3)  Starboard  Current 
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4)  Port  Voltage 

5)  Port  Current 


•  Starboard  AC  system  states 
6-8)  Voltages  a,  b ,  c 
9-11)  Generator  Currents  a,  b,  c 
12-14)  Propulsion  Currents  a,  b,  c 
15-17)  Power  Supply  Currents  a,  b,  c 

The  top  plot  in  Fig.  5-5  shows  the  normalized  integrated  variances  of  the  output 
states,  and  we  see  that  the  three  most  sensitive  outputs  to  the  triple  {1,2,3}  are 
the  power  supply  currents.  The  lower  plot  is  of  the  multi-dimensional  sensitivity 
(MS):  the  square  root  of  the  normalized  variance  divided  by  the  sum  of  the  variance 
of  the  uncertain  parameters  in  the  simulation.  This  tells  the  ratio  of  the  standard 
deviation  of  the  measured  states  to  the  standard  deviation  of  the  uncertain  inputs. 
There  is  an  eightfold  increase  in  the  variation  of  the  power  supply  currents,  which 


Output  parameter  number 
Multi-dimensional  sensitivity  parameter 


Output  parameter  number 


Figure  5-5:  Starboard  AC  system  normalized  integrated  variances  and  multi¬ 
dimensional  sensitivity  (MS)  for  ONR  IPS. 
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is  not  very  reassuring.  But  this  is  only  for  ±1%  uncertainty.  Fig.  5-6  shows  the 
integrated  variances  (above)  and  multi-dimensional  sensitivities  (below)  of  the  power 
supply  currents  for  ±1  —  17%  uncertainty.  The  near  eightfold  increase  only  occurs  at 
±1%  uncertainty,  and  stays  below  2  past  ±5%  uncertainty  (although  it  does  appear 
to  be  slightly  on  the  rise  past  ±15%). 


State  1 5 


State  1 6 


State  1 7 


%  individual  parameter  variation 


0.05  0.1  0.15 


0.05  0.1  0.15 


%  individual  parameter  variation 


Figure  5-6:  Multi-dimensional  sensitivity  (MS)  of  starboard  power  supply  currents 
for  ±1  —  17%  uncertainty  for  ONR  IPS. 


5.1.3  Adaptive  Simulation 

We  now  interpret  the  results  of  a  two  dimensional,  dimension  adaptive,  four  element 
simulation  of  the  uncertain  pair  {1,2},  the  bus  shunt  and  filter  capacitances.  Fig. 
5-7  shows  the  behavior  of  the  dimension  adaptive  algorithm  based  on  the  integral 
convergence  of  state  15,  the  a  power  supply  current. 

We  can  see  from  the  small  amount  of  computational  effort  expended  that  the 
model  is  a  relatively  well  behaved  function  of  the  bus  shunt  capacitance  below  its 
mean,  but  that  above  a  higher  resolution  is  needed  for  an  accurate  result  (Fig.  5-7). 
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Figure  5-7:  Record  of  dimension  adaptive  behavior  in  each  element  for  a  two  dimen¬ 
sional,  four  element  simulation  with  uncertainty  in  the  bus  shunt  and  filter  capac¬ 
itances  for  ONR  IPS.  Each  location  on  the  plot  represents  a  difference  grid  in  the 
corresponding  element.  Recently  evaluated  difference  grids  are  denoted  by  +  and 
old  difference  grids  by  O.  The  title  of  each  quadrant  is  the  portion  of  the  integral’s 
domain  covered  by  that  element.  It  also  corresponds  to  the  percent  variation  from 
the  mean  of  the  uncertain  couple. 


Perhaps  the  most  useful  observation  is  the  presence  of  higher  order  joint  difference 
grids  in  the  [0  1]  x  [—1  0]  element,  which  indicates  that  there  is  significant  coupled 
behavior  when  the  shunt  capacitance  and  the  filter  capacitance  are  respectively  above 
and  below  their  means.  This  is  not  a  measure  of  sensitivity;  however,  this  sort  of 
coupling  as  a  function  suggests  that  high  multidimensional  sensitivity  to  uncertainty 
in  a  set  of  parameters  may  be  a  result  of  the  interaction  between  them. 

Similar  information  can  be  inferred  from  higher  dimensional  simulations;  we  only 
show  a  two  dimensional  result  because  it  is  easily  presented  graphically.  The  same 
mean  and  variance  analysis  seen  in  the  previous  section  can  be  obtained  with  dimen¬ 
sion  adaptive  simulation,  but  we  omit  them  here  to  avoid  redundancy. 
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5.2  Pulse  Power  System 


We  analyze  uncertainty  in  a  Simulink  model  describing  the  operation  of  a  large  pulse 
load  reflecting  the  power  consumption  of  a  rail  gun  [8] .  The  alternator  is  charged  by 
accelerating  its  rotor  to  18,000  rpm,  at  which  point  the  inverter  and  charging  motor 
are  disconnected  from  the  alternator  and  a  shot  is  fired. 

A  six  dimensional  dimension-adaptive  simulation  using  Clenshaw-Curtis  points 
was  run,  with  ±10%  uniform  uncertainty  in  the  alternator  held  winding  resistance, 
primary  generator  one  held  winding  leakage  inductance,  inductance  of  charging  motor 
switches,  and  charging  motor  phase  resistance,  and  ±5%  uniform  uncertainty  in  the 
”on”  configuration  output  rectiher  diode  resistance  and  charging  motor  excitation 
hux. 

Fig.  5-8  shows  mean  trajectories  with  standard  deviations  for  a  few  states,  and 
5-9  shows  histograms  at  t  =  30  seconds.  A  total  of  53  system  evaluations  were 
required  to  attain  an  error  tolerance  of  10-3  for  the  charging  motor  phase  one  stator 
current.  In  most  dimensions,  the  algorithm  expended  little  effort,  the  highest  order 
joint  grid  evaluated  being  level  four  in  the  dimension  corresponding  to  the  charging 
motor  excitation  hux  and  level  three  in  that  corresponding  to  the  diode  resistance, 
indicating  some  coupling  between  those  parameters.  A  level  three  full  grid  would  have 
used  729  points,  taking  much  more  time  to  achieve  a  comparable  level  of  accuracy. 

It  can  be  seen  from  Fig.  5-8  that  the  generator  held  winding  voltage  and  generator 
speed  are  sensitive  to  the  simulated  uncertainty,  the  hrst  during  charging  before  the 
shot  is  hred  and  the  latter  afterward.  The  distributions  at  30  seconds  appear  to  mostly 
be  Gaussian,  with  the  exception  of  the  generator  speed,  which  is  nearly  uniformly. 


5.3  RTDS 

Uncertainty  analysis  was  performed  on  a  pulse  load  charging  event  for  a  notional  all- 
electric  ship  at  the  Center  for  Advanced  Power  Systems  at  Florida  State  University 
using  the  dimension-adaptive  algorithm  (section  3.1.2,  [13]).  The  system,  a  schematic 
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Charging  motor  phase  1  stator  current 


Primary  generator  voltage 


Figure  5-8:  Mean  trajectories  with  standard  deviation  envelopes  of  pulse  power  sys¬ 
tem  states. 
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Charging  motor  phase  1  stator  current 
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Figure  5-9:  Histogram  of  pulse  power  system  states  at  t  =  30  seconds. 
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of  which  is  given  in  Fig.  A-l,  has  two  36  MW  main  generators  and  two  4  MW  auxiliary 
generators.  The  load  is  comprised  of  two  36.5  MW  propulsion  motors,  and  maximum 
4  MW  ship  service  loads.  A  thorough  description  can  be  found  in  [25,40].  The 
notional  system  simulation  was  run  on  a  real  time  digital  simulator  [24],  on  which 
large  power  system  simulations  can  be  run  in  real-time  with  order  50 fxs  time  steps, 
and  into  which  physical  hardware  can  be  incorporated.  The  simulations  carried  out 
for  this  analysis  were  run  slower  than  real-time  with  a  15/rs  time  step. 

Fig.  A-2  in  appendix  A  shows  interpolated  versus  actual  values  for  the  maximum 
deviation  from  a  nominal  frequency  for  main  generator  one  and  auxiliary  generator 
two.  Detailed  results  can  be  found  in  [26]. 


5.4  Discussion 

Collocation  was  used  to  assess  the  behavior  of  three  systems  under  mid- dimensional 
uncertainty.  There  were  significant  computational  savings  over  full  and  sparse  grids 
for  simulations  in  which  the  dimension-adaptive  algorithm  was  used,  as  was  evidenced 
by  the  uneven  sampling  of  those  systems’  random  spaces.  Instructions  for  using  a 
Matlab  implementation  of  the  collocation  algorithms  is  given  in  appendix  B. 
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Chapter  6 


Unscented  Kalman  Filtering  with 
Dimension- Adaptive  Numerical 
Integration 


In  this  chapter,  the  dimension-adaptive  algorithm  is  implemented  in  an  unscented 
Kalman  filter.  The  performance  of  the  new  filter  is  evaluated  on  three  example 
systems  from  the  literature. 


6.1  Background 

We  want  to  estimate  the  state  xk  €  9ft”  of  the  nonlinear  discrete  dynamical  system 
with  process  and  measurement  noise 


Xfc+i  =  f(xfc,  uk,  k )  +  wk 

(6.1) 

zk  =  h(xfc,  Ufc,  k)  +  rk, 

where  w  and  r  are  zero-mean,  uncorrelated  noise  vectors  with  respective  covariances 
Q  and  R. 
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6.1.1  The  Extended  Kalman  Filter 


The  Kalman  filter  produces  optimal  state  estimates  based  on  current  and  previous 
measurements  by  first  computing  a  prediction  xfc+1|fc  via  a  system  model  and  then  a 
corrected  prediction  Xfc+i|fc+i  with  a  measurement  from  the  actual  system.  Both  the 
mean  and  covariance  matrix  P  of  the  state  are  propagated  with  the  filter.  When  f 
and  h  in  (6.1)  are  linear  with  system  matrices  A.  B,  C,  and  D,  the  predictor  portion 
is  calculated  by 

Xfc-fi|fc  Ax^i^  T  Bufc 

Pk+i\k  =  APfc|fcAT  +  Q  (6.2) 

^fc+i  Cx^_|_x|/c  Du^ 


with  correction 

Pot  =  CPfc+i|fcCT,  Pxz  =  Pfc+i|fcCT 

L  =  P xz l  (P vv  +  K) 

■^fe+i|fc+i  T  L(zfc+1  z^,_(_i) 

Pfc  +  l|fc+l  =  Pfc+l|fc  l^IZ' 


(6.3) 


In  the  nonlinear  case,  the  mean  and  output  can  be  calculated  in  the  same  fashion 
with  the  exact  system  model,  but  there  is  no  analogous  method  for  the  covariance 
matrix,  so  an  approximation  must  be  used  instead.  Traditionally,  this  has  been  done 
by  linearizing  about  x^ifc  and  then  performing  (6.2)  and  (6.3)  with  the  system  and 
output  Jacobians  F  and  H.  This  is  a  shortcoming  of  the  extended  Kalman  Filter: 
important  characteristics  of  the  model  are  lost  in  the  linearization,  and  the  Jacobian 
must  be  reevaluated  at  every  time  step,  a  tedious  and  computationally  expensive 
procedure. 


6.1.2  The  Unscented  Kalman  Filter 

In  [19],  the  necessity  for  linearization  in  the  extended  Kalman  filter  was  circumvented 
by  numerically  integrating  for  Pjufc,  Pzz,  and  Pxz.  The  original  system  equations  are 
evaluated  at  sigma  points  having  the  same  mean  and  covariance  as  the  state,  and  the 
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predicted  mean  and  covariance  are  then  calculated  from  the  transformed  points: 


cr  1  =  Xfc|fc,  wi  =  k/  (r)  +  hi) 

Oi  =  Xfcjfc  +  xk\k(y/(r)  +  K)Pfc|fc))«  Wi  =  1/2(77  +  k)  (6.4) 

<T+n  =  xfc|fc  -  Xfc|fc(^/ (77  +  n)Pk\k))i,  wi+v  =  1/2(77  +  «), 

where  7  =  2, . . . ,  77  +  1.  (A)j  denotes  the  ith  row  of  A.  and  the  matrix  B  with  square 
root  A  is  of  the  form  A7  A.  The  mean  and  covariances  are  then  calculated  by 

Xfc+i|fc  =  ES!1  Wif(<Ti),  zk+1  =  E ilV  7Uih(o-?:)  (6.5) 

Pk+i\k  =  E^l1  -  Xfc+i|fc)  (f(<Jj)  -  Xfc+i|jt)T 


Pxz  =  Ei=tX  «7*(f(t7i)  -  Xfc+l|fe)(h((Ti)  -  Zfc+1)T 


(6.6) 


P««  =  ESl1  W*(h(cr*)  -  Zfc_|_i)(h((7j)  -  Zfc+i)T. 


6.2  Dimension- Adaptive  Unscented  Kalman  Fil¬ 
tering 


The  sigma  points  of  [19]  are  a  type  numerical  integration,  specifically  one  that  assumes 
a  Gaussian  distribution,  as  with  Hermite  polynomial  Gaussian  quadrature.  It  is  clear 
that  the  any  valid  numerical  integration  scheme  could  be  used  in  place  of  the  sigma 
points,  so  it  is  sensible  to  employ  a  sophisticated  technique.  This  has  already  been 
done  to  some  extent:  [18]  implemented  the  unscented  filter  with  Gaussian  quadrature, 
and  [16]  with  quasi-Monte  Carlo.  We  propose  a  filter  that  uses  the  dimension-adaptive 
integration  scheme  from  [13]  in  place  of  the  sigma  points  of  the  original  unscented 
filter. 
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6.3  Performance  of  the  New  Filter 


We  now  compare  unscented  Kalman  filters  with  Kronrod-Patterson  dimension-adaptivity 
(DA),  Hermite  polynomial  Gaussian  quadrature  (GQ),  and  sigma  points  (SP)  and 
the  extended  Kalman  filter  (E)  on  the  falling  body  tracking  problem  from  [1],  the 
Kraichnan-Orszag  system  [30,42],  and  a  system  based  on  the  Lorenz  system  ex¬ 
ample  in  [18].  Gaussian  quadrature  with  Hermite  polynomials  integrates  functions 
over  a  Gaussian  distribution,  and  is  slightly  superior  to  Legendre  polynomial  Gaus¬ 
sian  quadrature  for  this  application.  However,  there  is  no  standard  nested  Gaus¬ 
sian  distributed  quadrature  sequence,  which  is  why  we  use  the  uniformly  distributed 
Kronrod-Patterson  quadrature  in  the  dimension  adaptive  filter. 

Most  of  the  previous  literature  present  results  via  the  absolute  value  of  the  mean 
error  for  all  trials  alongside  the  error  variance.  For  concision  and  because  the  relevant 
information  is  retained,  we  alternatively  give  the  mean  of  the  absolute  value  of  the 
error  in  our  experiments. 

In  each  experiment  100  Monte  Carlo  simulations  were  performed.  In  this  sections 
figures,  EKF  stands  for  the  extended  Kalman  filter,  LIKF  the  unscented  Kalman  filter, 
GQKF  the  Gaussian  quadrature  Kalman  filter  of  [18],  and  DAQKF  the  dimension- 
adaptive  quadrature  Kalman  filter. 


6.3.1  Falling  Body  Problem 

The  equations  of  motion  for  the  system  are  given  by 


Xi  =  —x2  +  w  i 


x2  =  -  exp(— 7£i)x|x3  +  w2  (6.7) 

x3  =  w3, 


with  measurement 


z  =  \J  M2  +  {x\  —  H)2  +  r, 


(6.8) 


where  7  =  5-  10~5,  M  =  105,  H  =  105,  and  w  and  r  are  zero-mean  uncorrelated 
noises  with  covariances  given  by  Q  =  0  and  R  =  104,  respectively.  The  process  noise 
is  removed  for  the  reason  that  it  can  mitigate  linearization  errors  and  thus  actually 
improve  the  the  extended  Kalman  filter’s  performance.  The  initial  state  of  the  true 
system  is 

x(0)  =  [  3  .  io5  2  •  104  10~3  ]T  (6.9) 

The  initial  estimates  for  the  simulated  filter  were 

x(0)  =  [  3 • 105  2 • 104  3 • 10"5  }T 

106  0  0  (6.10) 

P0!o  =  0  4  •  106  0 

0  0  10”4 

With  the  faulty  initial  estimate  we  are  assuming  the  body  to  be  heavier  than  it 
actually  is.  Measurements  are  made  at  every  second,  with  64  fourth-order  Runge- 
Kutta  system  updates  between  each  measurement. 

Little  to  no  improvement  is  seen  over  the  Gaussian  quadrature  and  sigma  point 
Liters  (Fig.  6-1). 


6.3.2  The  Kraichnan-Orszag  System 

The  Kraichnan-Orszag  three  mode  system  can  be  highly  oscillatory,  depending  on  ini¬ 
tial  conditions,  and,  with  infrequent  measurements,  represents  a  challenging  filtering 
and  stochastic  simulation  problem.  The  system  equations  are 

Xi  =  x2x3  +  w  1 

X2=X1X3+W2  (6-H) 

x3  =  —2xix2  +  w3. 
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Figure  6-1:  Mean  absolute  error  for  position,  velocity,  and  £3  for  falling  body  problem. 


As  with  the  falling  body  problem,  the  process  noise  covariance  matrix  is  set  Q  =  0. 
For  the  measurement  we  use 

z  =  x  +  r,  (6-12) 

with  measurement  noise  covariance  R  =  51.  The  initial  state  used  in  both  the  ’true’ 
and  filter  simulation  of  (6.11)  is  x(0)  =  [  0  1  20  ]T,  with  initial  covariance  estimate 
P0|o  =  I.  A  fourth  order  Runge-Kutta  solver  with  dt  =  0.01  was  used,  with  measure¬ 
ments  taken  every  every  20 dt.  Fig.  6-2  shows  the  corresponding  state  trajectories. 
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Three  dimensional  trajectory 


States  1 , 2,  3  versus  time 


Figure  6-2:  Trajectories  of  the  Kraichnan-Orszag  system. 

Fig.  6-3  shows  the  errors  of  dimension- adaptive  filters  with  different  error  toler¬ 
ances.  The  filters  with  smaller  tolerances  outperform  those  with  larger  tolerances. 
Initially,  more  points  are  used  by  the  adaptive  algorithm  to  compute  more  precise 
means  and  covariances.  As  the  filters’  confidences  grow,  the  covariances  shrink,  de¬ 
creasing  the  size  of  the  region  being  integrated  over.  Numerical  integrals  are  easier 
on  smaller  domains,  and  hence  fewer  points  are  used  by  the  adaptive  filter  as  the 
covariance  decreases. 

These  errors  are  very  choppy;  the  reason  for  this  is  in  the  nature  of  the  Kraichnan- 
Orszag  system.  Fig.  6-4  shows  a  ’true’  trajectory  in  black,  and  noisy  trajectories 
with  noise  added.  The  noisy  trajectories  sometimes  take  wrong  turns  at  the  corners, 
causing  the  error  to  spike,  and  then  return  to  the  correct  path  at  the  next  corner, 
causing  the  error  to  spike.  Once  the  variance  has  decreased  sufficiently,  the  points 
are  close  enough  together  that  they  all  take  the  same  path,  and  the  error  smoothens 
out. 

The  dimension-adaptive  filter  outperforms  all  of  the  other  filters  on  this  exam¬ 
ple  (fig.  6-5).  Initially,  it  expends  a  significantly  larger  effort  than  the  others,  but 
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State  1  mean  absolute  error 


State  2  mean  absolute  error 


Time 


Mean  number  of  points  propagated 


0  5  10 

State  3  variance 


0  5  10 

Time 


Figure  6-3:  Error,  variance,  and  mean  number  of  points  propagated  for  dimension- 
adaptive  filters  with  different  error  tolerances  on  Kraichnan-Orszag  system. 
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Figure  6-4:  Noisy  Kraichnan-Orszag  trajectories. 


73 


eventually  only  requires  about  the  same  effort  as  the  unscented  filter  to  achieve  its 
prescribed  error. 


Time 


Figure  6-5:  Performance  of  filters  on  the  Kraichnan-Orszag  system. 


6.3.3  The  Lorenz  system 

Many  systems  of  interest  follow  trajectories  like  those  of  the  Lorenz  system,  e.g. 
electrical  circuits  [5]  and  lasers  [44],  A  comparison  of  is  shown  on  an  example  based 
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on  the  Lorenz  system  from  [18].  The  system  equations  are  given  by 


X\  =  cr(— Xi  +  x2) 

x2  =  px  1  —x2  —  xix 3  (6.13) 

x3  =  ~/3x  3  +  XiX2  +  0.5w, 

where  a  =  10,  p  =  28,  and  [3  =  8/3,  which  correspond  to  the  well-known  butter¬ 
fly  effect  system.  Fig.  6-6  show  system  trajectories  for  the  parameters  and  initial 
conditions  used.  The  output  is  given  by 

States  1,2,3  versus  time  Three  dimensional  trajectory 


Figure  6-6:  Trajectories  of  the  Lorenz  system. 


y  =  yj(xi  -  0.5)2  +  x\  +  x\  +  0.065r,  (6.14) 

and  the  noises  w  and  r  each  have  the  same  variance  Q  =  R  =  dt. 

An  Euler  step  with  dt  =  0.01  is  used.  In  [18],  measurements  were  taken  at  every 
time  step;  here  instead  we  take  measurements  at  10 dt. 

The  dimension-adaptive  filter  outperforms  the  other  unscented  Liters,  avoiding 
the  first  three  error  spikes  (Fig.  6-7).  The  number  of  points  used  goes  above  and 
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Error  Error 


below  the  27  used  by  the  Gauss-Hermite  filter,  depending  on  the  number  of  points 
needed  to  achieve  accurate  integration  of  the  state  and  covariance.  The  extended 
Kalman  filter  fails  on  this  example,  and  so  it  is  not  shown. 


State  1  mean  absolute  error 


15 
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State  2  mean  absolute  error 
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State  3  mean  absolute  error 
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Figure  6-7:  Performance  of  filters  on  the  Lorenz  system. 
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6.4  Summary 


A  new  filter  using  dimension-adaptive  numerical  integration  within  a 
unscented  Kalman  filter  is  presented.  Performance  gains  were  seen  on 
systems. 


conventional 
two  example 
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Chapter  7 


Summary 


Stochastic  differential  equations  provide  a  framework  for  modeling  systems  with  un¬ 
certainty  and  informing  design  decisions,  but  for  large-scale  systems  they  are  of¬ 
ten  hard  to  solve  numerically  because  of  high  dimensional  random  spaces.  Collo¬ 
cation  is  an  attractive  approach  for  large-scale  systems  with  moderate  dimensional 
ranges  because  of  its  ’black-box’  implementation  and  the  existence  of  sophisticated 
tensor-product  algorithms  with  which  it  can  be  used.  In  this  work,  full,  sparse,  and 
dimension-adaptive  [13]  grid  collocation  with  various  polynomial-based  component 
quadratures  were  compared  on  standard  test  integrals  [11]  and  then  applied  to  com¬ 
pute  solutions  to  stochastic  differential  equations  describing  a  notional  electric  ship 
in  standard  operation,  a  pulse  power  system,  and  a  charging  event  on  a  notional 
electric  ship.  Statistical  moments  and  PDF’s  were  extracted  from  these  solutions, 
from  which  sensitivity  was  inferred.  The  capability  to  assess  the  behavior  of  a  system 
under  uncertainty  is  a  powerful  design  tool.  It  is  an  indicator  of  robust  design  and 
control  strategies;  in  the  case  of  the  electric  ship,  two  such  examples  are  safe  load¬ 
ing  conditions  or  ranges  of  deviation  from  nominal  behavior  as  a  result  of  imprecise 
system  knowledge.  A  Matlab  toolbox  has  been  created  and  released  to  the  Electric 
Ship  Research  and  Development  Consortium  Partners  (appendix  B).  The  results  of 
Chapter  5  were  generated  with  this  toolbox;  there  are  also  additional  users  at  this 
writing. 

The  sparse  grid  with  shifted  indexing,  a  modification  to  the  conventional  sparse 
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grid  formulation,  was  presented,  along  with  justification  in  the  form  of  theoretical 
bounds  and  numerical  experiments.  The  dimension-adaptive  grid  collocation  algo¬ 
rithm  was  then  applied  to  unscented  Kalman  filtering,  and  improvement  over  the 
extended  and  other  unscented  Kalman  filters  was  seen  for  two  nonlinear  filtering 
problems.  An  adaptive  unscented  filter  is  promising  because  of  its  potential  for  com¬ 
putational  savings;  in  general  less  computational  effort  is  required  as  the  filter  attains 
greater  confidence  in  its  estimates.  The  gains  for  the  examples  shown,  however,  are 
modest. 

Future  work  includes  collocation  using  locally  adaptive  grids,  i.e.  algorithms  that 
can  adapt  to  regionally  rather  than  dimensionally  in  the  random  space.  Another 
direction  is  the  fusion  of  Monte  Carlo  and  collocation.  While  at  first  seemingly  counter 
to  what  makes  Monte  Carlo  work,  which  in  part  is  that  points  are  placed  directly 
into  the  multi-dimensional  random  space  rather  than  built  up  from  one-dimensional 
bases,  assigning  certain  ’easier’  dimensions  to  collocation  and  the  rest  to  an  aggregate 
Monte  Carlo  dimension  may  yield  faster  convergence,  and  be  particularly  effective  if 
done  adaptively. 
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Appendix  A 


Figures  from  FSU  CAPS 
Collaboration 


Zone  4 
—  I  PCM4 


—  /I  Zone  2 
PCM4 


Prop 
Drive  1 


Pulse 

Load 


G)  (J3; 

MTG1  ATG1 


Figure  A-l:  Schematic  for  RTDS  notional  all-electric  ship  at  Florida  State  University. 
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Maximum  Frequency  Deviation  (pu) 


Value  Predicted  by 
Surrogate  Model 


o 

o 


O  MTG1 
O  ATG2 


Actual  Value 


0.00 


Figure  A-2:  Comparison  between  interpolated  and  true  values  for  maximum  frequency 
deviation  of  main  generator  one  and  auxiliary  generator  two  for  RTDS  notional  elec¬ 
tric  ship  at  Florida  State  University. 
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Appendix  B 


Readme  for  Matlab 
implementation  of  Collocation 
Algorithms 

B.l  Overview  and  introduction 

This  appendix  serves  as  a  guide  for  the  Matlab  program  Collocation. m ,  a  set  of 
routines  for  numerically  integrating  and  interpolating  multi-dimensional  functions  on 
the  unit  cube  via  tensor  products  of  one  dimensional  rules.  It  can  be  downloaded 
from  http://web.mit.edu/hovergroup/docs.html.  Monte  Carlo  and  quasi-Monte  Carlo 
techniques  are  not  part  of  its  functionality,  although  in  many  situations,  namely  when 
dealing  with  high-dimensional  or  low-smoothness  functions,  they  are  the  best  choice. 
The  approaches  included  are: 

•  Full  grid  tensor-product  cubature 

•  Sparse  grid  tensor-product  cubature  [12,28,38] 

•  Uniform-adaptive  full  or  sparse  grid  tensor-product  cubature 

•  Dimension-adaptive  tensor-product  quadrature  [13] 

•  Uniform  subspace  splitting 
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Each  approach  is  based  on  one- dimensional  quadratures.  Polynomial  based  rules  have 
the  highest  polynomial  exactness,  which  is  the  highest  degree  polynomial  a  method 
will  integrate  exactly  using  a  fixed  number  of  points.  The  available  one-dimensional 
basis  quadratures  are: 

•  Legendre  polynomial  Gaussian  quadrature  [36] 

2 n  —  1  polynomial  exactness,  not  nested,  constant  integration  kernel. 

•  Hermite  polynomial  Gaussian  quadrature 

2 n  —  1  polynomial  exactness,  not  nested,  Gaussian  integration  kernel. 

•  Clenshaw- Curtis  [4,41] 

n  polynomial  exactness,  nested,  constant  integration  kernel. 

•  Kronrod-Patterson  [10,31,33] 

2 n  +  m  —  1  polynomial  exactness,  nested,  constant  integration  kernel.  Kronrod- 
Patterson  is  a  nested  sequence  of  quadratures  beginning  with  Gaussian  quadra¬ 
ture.  n  is  the  number  of  Gaussian  quadrature  points,  and  m  the  number  of 
extension  points.  Specifically,  the  sequence  beginning  with  the  three  points 
Legendre  polynomial  rule  is  included.  The  Kronrod-Patterson  rule  is  not  gen¬ 
erated  by  a  program,  but  read  from  a  table  that  goes  up  to  the  seventh  level 
in  a  sequence.  For  this  reason,  it  is  important  when  using  an  adaptive  grid  to 
set  the  maximum  grid  level  below  the  highest  available  Kronrod-Patterson  rule; 
otherwise  a  higher  rule  will  be  called,  and  the  program  will  return  an  error. 

Gaussian  quadrature  has  the  highest  polynomial  exactness,  but  it  is  not  nested,  and 
for  this  reason  Clenshaw- Curtis  and  Kronrod-Patterson  are  superior  choices  for  sparse 
grids  or  adaptive  integration. 

B.2  Using  Collocation.m 

Collocation :m  is  called  by 

[outputs]  =  Collocation(q ,  d,  @ integrand ,  options ) 
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where  q  is  the  level  and  d  the  dimension,  q  is  only  used  as  the  non-adaptive  default 
level  if  no  adaptive  options  are  specified.  @ integrand  is  the  function  handle  for  the 
integrand  function.  The  outputs  correspond  to  the  solution  to  the  integral  plus  what 
is  specified  in  the  options,  for  example,  if  ../ numval 1...  is  the  only  input  after  inte¬ 
grand,  the  output  will  be  [y,n],  where  y  is  the  computed  value  of  the  integral  and  n 
the  number  of  function  evaluations. 

options : 

•  grid, 

~  CC:  Clenshaw-Curtis 

—  KPL3:  Kronrod-Patterson 

—  Leg:  Legendre  polynomial  Gaussian  quadrature 

—  ffer:  Hermite  polynomial  Gaussian  quadrature 

•  s-f 

—  0:  full  grid 

—  1:  sparse  grid 

•  funin 

—  Additional  inputs  to  the  integrand  function  in  addition  to  the  evaluation 
point. 

•  timer 

—  1:  Output  the  time  taken. 

•  showgrid 

—  1:  In  2  dimensions,  display  the  points  and  graph  of  the  function  as  they  are 
evaluated.  In  the  dimension  adaptive  case,  the  old  and  active  index  sets 
for  the  current  element  are  also  shown.  In  3  dimensions,  the  evaluation 
points  are  plotted. 
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•  variance 


—  1:  Calculate  the  variance. 

•  mea  (Multi- Element  Adaptivity) 

—  Inputs  are  specified  as  an  array: 

1.  Number  of  partitions  (1  if  none) 

2.  Desired  resolution  for  each  grid  element.  For  non-adaptive,  input  an 
integer.  That  number  will  be  the  level  in  each  partition. 

3.  Maximum  grid  level  (to  prevent  the  adaptivity  from  getting  out  of 
control) 

4.  Adaptivity  -  0:  non  or  uniform  adaptivity  (depending  on  2),  1:  Di¬ 
mension  adaptivity 

5.  1:  Return  record  of  adaptive  behavior 

6.  Adaptive  criterion  for  time  dependence:  0:  Integral,  1:  Maximum 

7.  For  integrands  with  vector  ranges,  this  specifies  the  output  parameter 
the  adaptivity  responds  to  (only  for  time  dependent). 

•  numval 

—  1:  Return  the  number  of  function  evaluations. 

•  timedep 

—  If  the  integrand  returns  vectors  in  time,  this  input  should  be  a  two  element 
array  in  the  form  [dt,tf\,  where  dt  is  the  time  step  the  simulation  will  be 
splined  to  and  tf  the  end  time  of  the  simulation. 


•  nout 


—  In  the  time  varying  case,  this  is  the  number  of  function  outputs. 


•  interp 


—  An  n  by  d  matrix,  where  n  is  the  number  of  points  being  interpolated,  d  the 
dimension,  and  each  row  a  point  being  interpolated.  Do  not  use  subspace 
splitting  or  Hermite  polynomial  Gaussian  quadrature  if  using  this  feature. 

•  show  prog 

—  1:  Display  which  grid  is  being  evaluated  and  the  current  attained  error 
tolerance  (dimension-adaptivity  only). 

The  function  /  should  be  of  the  form 

function  y  =  f(x,  funin )  for  scalar  functions. 

function  [y,t]  =  f(x,  funin )  for  functions  that  are  vectors  through  time. 

t  is  the  vector  of  times  the  solutions  is 
splined  to,  of  the  form  t\  :  dt  :  £2- 

where  a:  is  a  vector  of  the  coordinates  for  the  point  to  be  evaluated  and  funin  any 
extra  information  for  the  function.  If  there  are  no  additional  inputs  to  the  function, 
keep  an  unused  second  input.  It  should  return  the  function  value  at  that  point, 
which  can  be  a  scalar  or  a  vector(s)  through  time.  A  description  of  the  output  of  / 
is  provided  to  the  program  through  the  options. 

Sometimes  there  are  multiple  values  of  interest;  this  is  often  the  case  in  large 
scale  dynamic  simulation,  when  there  are  many  states.  Because  each  function  eval¬ 
uation  can  be  costly,  it  is  sensible  to  be  able  to  integrate  each  state  using  one  set  of 
simulations.  In  this  case  multiple  states  can  be  evaluated,  and  a  matrix  with  rows 
corresponding  to  states  and  columns  to  times  is  returned. 

When  function  evaluations  are  very  expensive,  it  is  highly  recommended  to  have 
Collocation.m  call  an  outer  function  which  evaluates  and  writes  to  a  hie  the  function 
value  if  it  is  a  new  evaluation  point,  and  reads  the  function  value  from  an  existing 
hie  if  it  has  already  been  evaluated  at  that  point.  It  may  happen  that  a  simulation 
is  stopped  before  all  the  grids  have  been  computed,  or  more  similar  simulations  are 
run  in  the  future;  writing  the  evaluations  to  a  hie  as  the  program  goes  prevents  lost 
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or  repetitive  calculation.  Below  is  an  example  of  how  such  a  program  might  look. 


function  y  =  outer  (x,  funin ) 

if  /  has  previously  been  evaluated  at  x 
read  f(x,  funin )  from  hie 

else 

evaluate  f(x,  funin ) 
write  f(x ,  funin)  to  a  hie 
end 

V  =  f(x,  funin) 

B.2.1  A  warning  about  computational  costs 

Don’t  set  the  adaptive  tolerance  too  small  or  the  level  of  the  non-adaptive  grid  too 
high.  Computational  cost  can  increase  exponentially:  if  a  certain  level  takes  an  hour, 
the  next  level  up  could  take  a  year.  For  the  adaptive  grids,  there  is  an  option  to  impose 
a  maximum  grid  level.  It’s  easy  to  get  a  feel  for  how  long  the  program  takes  with  the 
test  integrals  (which  are  more  difficult  to  accurately  integrate  than  the  ODE). 


B.3  Description  of  each  program 

•  Collocation  hies 

—  C  allocation. m 

Interface  program  for  user.  Assembles  options  (see  below)  and  calls  sub¬ 
programs. 

—  ColOutput.m 

Calls  multi-dimensional  integration  programs  and  assembles  computed  re¬ 
sults. 

—  C  olP  artition.m 

Partitions  the  region  of  integration  into  identical  subspaces. 


ColFG.m 


Generates  and  evaluates  a  full  grid  tensor-product  quadrature. 

ColSG.m 

Generates  and  evaluates  a  sparse  grid  tensor-product  quadrature. 

ColSG  Index  Sets. m 

Creates  the  index  sets  for  the  sparse  grid. 

ColDif  fGrid.m 

Returns  abscissas  and  weights  for  one-dimensional  quadrature  routines.  If 
specified,  returns  one  dimensional  difference  grid. 

ColS  ortComb .  m 

Combines  coincident  points  in  sparse  grids. 

ColAdaptUnif  :m 

Uniform  adaptivity:  evaluates  higher  order  grids  until  adaptive  tolerance 
is  reached. 

Col  Adapt  Dim. m 

Dimension  adaptivity:  the  dimension  adaptive  algorithm  from  [13]. 

Col  Adapt  DimGr  id. m 

Generates  the  grid  for  Col  Adapt  Dim.m. 

Col  Adapt  DimV  alid.m 

Checks  that  the  next  grid  evaluated  in  ColAdaptDim.m  is  valid. 
ColInterpFG.m 

Interpolates  chosen  interpolation  points  if  a  full  grid  has  been  used. 

C olInterpSGD  A.m 

Interpolates  chosen  interpolation  points  if  a  sparse  of  dimension-adaptive 
grid  has  been  used. 

ColInterpDif  fb.m 

Generates  the  one-dimensional  difference  coefficients  for  constructing  sparse 
and  dimension-adaptive  grid  interpolants. 


—  FindLegendreRoots.m  (from  Numerical  Recipes  [36]) 

Returns  the  abscissas  and  weights  for  Legendre  polynomial  Gaussian  quadra¬ 
ture. 

—  HermQuad.m 

Returns  the  abscissas  and  weights  for  Hermite  polynomial  Gaussian  quadra¬ 
ture. 

—  ClenshawCurtis.m 

Returns  abscissas  and  weights  for  Clenshaw- Curtis  quadrature. 

—  KPL3:m 

Returns  abscissas  and  weights  for  Gauss-Kronrod- Patterson  quadrature 
beginning  with  the  three  point  Legendre  polynomial  Gaussian  quadrature. 


•  Auxiliary  hies 


—  Testlntegrand.m 

Test  functions,  including  those  in  [11], 

—  TestlntegrandSol.m 

Exact  solutions  to  test  functions  in  Testlntegrand.m. 

—  Plots JD2.m 

Plots  two-dimensional  instantiations  of  each  test  function. 

—  dirinit.m 

Adds  the  directory  to  the  Matlab  path,  so  Collocation,  can  be  called  from 
any  location. 

—  CollocationDemo.m 

Runs  Collocation.m  on  test  functions  and  a  first  order  ODE. 
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B.4  Test  functions 


The  ten  test  functions  in  Testlntegrals.m  are  shown  below.  The  first  six  are  from 

[11]- 


OSCILLATORY 

/(x) 

=  cos(27TUi  +  Eil  CiiXi) 

CORNER  PEAK 

/(x) 

=  (1  +  Eti  aiXi)-^ 

PRODUCT  PEAK 

/(x) 

=  n:L1(a-2  +  (^-ui)2)-1 

GAUSSIAN 

/(x) 

=  exp(—  Eli  a2i(xi  -  Ui)2) 

C° 

/(x) 

=  exp(—  Eti  Oil^i  -  Ui\) 

DISCONTINUOUS 


UNCOUPLED 


MAX  PRODUCT 


ODD  C° 


/(x)  = 


0  Xi  >  Ui  or  x2  >  u2 

exP(Ef=i  &iXi)  otherwise 


/(x)  =  Ef=i  exp  (difa  -  u^) 
/(x)  =  nti  ma x(x*  -Ui  +  1, 1) 
/(x)  =  exp  (-  Ez  odd^l^  - 


exp 


(  E i  even(a*(‘c*  Ui 


ODD  DISCONTINUOUS  /(x)  = 


0  3 i  odd  s.t.  Xi  >  Ui 

exp(Ef=i  ciiXi)  otherwise 

*a  and  u  are  tunable  parameters. 


All  entries  of  a  are  positive  and  u  is  randomly  generated.  Run  the  program 
PlotS-D2:m  to  see  the  two  dimensional  instantiation  of  each  function. 
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B.5  Running  CollocationDemo.m 


The  program  CollocationDemo.m  runs  the  package  on  some  examples.  Set  type  =  0 
to  integrate  and  interpolate  the  test  functions  above.  Set  type  =  1  to  evaluate  the 
expected  value  of  the  solution  to  the  ODE. 

X  =  ft  i  X  x(0)  =  /i'2 , 

where  [i\  and  /i2  are  uniformly  distributed  random  variables  on  the  interval  [0,1]. 
The  solution  to  the  ODE  is  x  =  This  is  treated  as  a  vector  valued  function 

in  which  each  entry  is  a  point  in  time.  Notice  that  the  expected  value  (what  is  being 
calculated) 

x  =  —(l-e~t) 

2C  ’ 

is  not  a  solution  to  the  ODE.  You’ll  see  that  the  dimension  adaptive  algorithm  requires 
a  much  smaller  adaptive  tolerance  to  achieve  a  good  error  for  the  variance;  this  is 
because  the  integrand  being  adapted  to  is  the  mean  and  not  the  variance. 

We  evaluate  the  electric  ship  in  the  same  fashion,  except  that  rather  than  having 
a  closed  form  solution  to  evaluate,  a  simulation  must  be  run  at  each  point  in  the 
integral’s  domain. 

For  sensitivity  analysis  we  look  at  the  variance  solution,  and  either  consider  the 
maximum  of  the  variance  through  time  or  the  integral  of  the  variance  through  time. 

Histograms  of  interpolated  points  are  generated  as  well;  the  usefulness  of  this 
feature  is  more  apparent  in  the  time  dependent  case,  where  straight  Monte  Carlo 
simulation  would  be  an  overly  expensive  means  of  finding  a  PDF  at  later  times. 
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